rm(list = ls())

library(foreign)
library(readstata13)
library(RcppRoll)
library(Rcpp)
library(tm)
library(quanteda)
library(list)
library(sandwich)
library(ggplot2)
library(dplyr)
library(hrbrthemes)


##########Load in Data for Studies 1, 2, and 3##########

#####Study 1#####
Qualtrics <- read.dta13("Blair_Chu_Schwartz_Study1_Data.dta")

#####Study 2#####
Turk <- read.dta13("Blair_Chu_Schwartz_Study2_Data.dta")

#####Study 3#####
Lucid <- read.dta13("Blair_Chu_Schwartz_Study3_Data.dta")


####################Main Text####################

##########Table 2##########

#####Study 1#####
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$chem2, na.rm=TRUE)
    B <- mean(temp$chem3, na.rm=TRUE)
    C <- mean(temp$chem_ask, na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Qualtrics)

##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])

#####Study 2#####
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==1)], temp$wtChem[which(temp$treatment==1 & temp$t_chem==1)], na.rm=TRUE)
    B <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==1)], temp$wtChem[which(temp$treatment==0 & temp$t_chem==1)], na.rm=TRUE)
    C <- weighted.mean(temp$direct[which(temp$t_chem==1)], temp$wtChem[which(temp$t_chem==1)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.3 <- bootTreat2(Turk)

##Indirect Support##
round(c(full.3$Indirect, quantile(full.3$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,1] > 0))/length(full.3$boot[,1])

##Direct Support##
round(c(full.3$Direct, quantile(full.3$boot[,2], c(0.025, 0.975))), digits=4)
#P-Value#
1-length(which(full.3$boot[,2] > 0))/length(full.3$boot[,2])

##Difference (Indirect-Direct)
round(c(full.3$Superficial, quantile(full.3$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,3] > 0))/length(full.3$boot[,3])


##########Results Mentioned In-Text##########

#####Mean Number of Policies Supported in the Control List#####

##Study 1##
summary(Qualtrics$chem3)

##Study 2##
weighted.mean(Turk$list[Turk$treatment==0 & Turk$t_chem==1], Turk$wtChem[Turk$treatment==0 & Turk$t_chem==1], na.rm=TRUE)

#####Testing the Difference in List-Direct Support Between Study 1 and 2 (P-Value)#####
1-length(which(full.2$boot[,3] < full.3$boot[,3]))/length(full.2$boot[,3])

#####Testing the Difference in Direct Support Between Study 1 and 2 (P-Value)#####
1-length(which(full.2$boot[,2] > full.3$boot[,2]))/length(full.2$boot[,2])


##########Table 3##########

#####Study 1#####
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$chem2, na.rm=TRUE)
    B <- mean(temp$chem3, na.rm=TRUE)
    C <- mean(temp$chem_ask, na.rm=TRUE)
    D <- 1
    bootResults[i,1] <- (D) - (A - B)
    bootResults[i,2] <- (A-B) - (C)
    bootResults[i,3] <- C
    drop(list())
  }
  return(list(model=label, boot=bootResults, Sincere=mean(bootResults[,1], na.rm=TRUE), 
              Insincere=mean(bootResults[,2], na.rm=TRUE), 
              Non=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Qualtrics)

##Sincere Norm-Holders##
round(c(full.2$Sincere, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Insincere Norm-Holders##
round(c(full.2$Insincere, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Non Norm-Holders
round(c(full.2$Non, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])

#####Study 2#####
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==1)], temp$wtChem[which(temp$treatment==1 & temp$t_chem==1)], na.rm=TRUE)
    B <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==1)], temp$wtChem[which(temp$treatment==0 & temp$t_chem==1)], na.rm=TRUE)
    C <- weighted.mean(temp$direct[which(temp$t_chem==1)], temp$wtChem[which(temp$t_chem==1)], na.rm=TRUE)
    D <- 1
    bootResults[i,1] <- (D) - (A - B)
    bootResults[i,2] <- (A-B) - (C)
    bootResults[i,3] <- C
    drop(list())
  }
  return(list(model=label, boot=bootResults, Sincere=mean(bootResults[,1], na.rm=TRUE), 
              Insincere=mean(bootResults[,2], na.rm=TRUE), 
              Non=mean(bootResults[,3], na.rm=TRUE)))
}

full.3 <- bootTreat2(Turk)

##Sincere Norm-Holders##
round(c(full.3$Sincere, quantile(full.3$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,1] > 0))/length(full.3$boot[,1])

##Insincere Norm-Holders##
round(c(full.3$Insincere, quantile(full.3$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,2] > 0))/length(full.3$boot[,2])

##Non Norm-Holders
round(c(full.3$Non, quantile(full.3$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,3] > 0))/length(full.3$boot[,3])


##########Percentage of Non Norm-Holders Willing to Sign a Public Petition##########
Turk$Petition <- NA 
Turk$Petition[Turk$petition=="No"] <- 0
Turk$Petition[Turk$petition=="Yes"] <- 1

summary(Turk$Petition[Turk$t_chem==1])


##########Percentage of Friends/Family React Negatively to Supporting Chemical Weapons Use##########

#####Read in Data#####
OpenEnded <- read.csv("Blair_Chu_Schwartz_Taboo.csv")

#####Estimate Percentage of Friends/Family That React Negatively  to Supporting Chemical Weapons Use#####
Open_Friends <- subset(OpenEnded, open_friends!="")

Open_Friends$Negative <- ifelse(grepl("horribly|cruel|disagree|opposed|oppose|against|(even.*entertain)|(not.*support)|(wouldn’t.*support)|cruel|scary|(not.*believe)|(not.*agree)|(not.*use)|frightened|terrified|crazy|unwilling|outraged|(wouldn’t.*agree)|(should.*not)|upset|shouldn’t|hesitate|ostracized|(no.*agree)|(say.*no)|appalled|mad|reject|negative|negatively|hate|(I.*don't.*think.*any.*of.*my.*friends.*or.*family.*would.*support.*this)|(not.*a.*good)|(bad.*idea)|(Don’t.*think)|ashamed|(none.*of.*them)|(not.*at.*all)|inhumane|stupid|innocent|(not.*accept)|(not.*want)|(not.*be.*willing)|devastating|(don’t.*see)|reject|upset|aghast|fear|horrified|(wont.*like)|(won’t.*like)|horrible|WTF|(not.*be.*in.*favor)|disgust|protest|anger|(say.*no)|scared|sad|(do.*not.*need)|poorly|vehemently|horrific|unacceptable|unfavorably|anger|(all.*hate)|(find.*that.*drastic)|shouldnt|disappointed|(very.*bad)|harmful|(no.*circumstances)|(too.*dangerous)|horrors|harm|inhumane|(wouldn't.*like)|(wouldn't.*agree)|(would.*not.*agree)|condemn|(foolish.*idea)|abhorrent|insane|unethical|(is.*wrong)|(big.*no)|(not.*like)|(absolutely.*not)|(never.*be.*used)|worried|(more.*harm.*than.*good)|(not.*be.*used.*in.*any.*instance)|(not.*be.*supportive)|(not.*be.*okay)|(not.*be.*in.*favor)|(not.*be.*happy)|(opposing.*the.*idea)|(shake.*their.*heads)|banned|(freak.*out)|detestable|leery|weary|angry|(highly.*concerned)|dismayed|(won't.*support)|shocking|(not.*right)|(more.*harm)|(nobody.*will.*support)|evil|disapprove|(zero.*support)|(not.*well)|(not.*in.*favor)|(not.*good)|(not.*a.*good)|(no.*one.*I.*know)|(neither.*friends)|(neither.*family)|negetive|(kill.*children)|shocked|(wouldn't.*approve)|disgusted|alarmed|(too.*controversial)|outrage|(shouldn't.*happen)|(not.*ethical)|horror|disturbed|(not.*do.*it)|mortified|(never.*use)|(not.*favor)|hypocritical|(don't.*think.*they.*would.*like)|(don't.*think.*they.*would.*support)|terrible|devastating|repulsed|outrage|(I.*can't.*think.*of.*anyone.*I.*know)|(hell.*naw)|harshly|barbaric|badly|(attacking.*that.*policy.*position)|pissed|immoral|whoa|(not.*approve)|disapproval|apposed|(don't.*think.*any.*of)|(don't.*believe.*any.*of)|(no.*one.*would.*really.*like)|(most.*would.*not)|furious|(end.*of.*humanity)|(don't.*think.*they.*would.*agree)|opposition|(don't.*believe.*anyone.*I.*know)|(not.*happy)|(don't.*think.*my.*would.*agree)|(end.*of.*the.*world)|astonished|(object.*to.*it)|(no.*way)|(don't.*think.*the.*majority.*agree)|(don't.*.*think.*any.*would.*be.*in.*favor)|shock|disagree|(not.*approve)|(no.*need.*to.*use)|(too.*risky)|(thats.*wrong)|bad", Open_Friends$open_friends, ignore.case=TRUE), 1, 0)

summary(Open_Friends$Negative[Open_Friends$t_chem==1])


##########Nuclear Weapons Results##########

#####Study 2 (mTurk)#####
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==0)], temp$wtNuke[which(temp$treatment==1 & temp$t_chem==0)], na.rm=TRUE)
    B <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==0)], temp$wtNuke[which(temp$treatment==0 & temp$t_chem==0)], na.rm=TRUE)
    C <- weighted.mean(temp$direct[which(temp$t_chem==0)], temp$wtNuke[which(temp$t_chem==0)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Turk)

##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])


#####Study 3 (Lucid)#####
##Number of List Items Chosen Among those in Either of the Two Lists 
Lucid$controllist <- 0
Lucid$controllist[Lucid$nonsensitivelist==1] <- 1
Lucid$controllist[Lucid$nonsensitivelist==2] <- 2
Lucid$controllist[Lucid$nonsensitivelist==3] <- 3

Lucid$treatlist <- 0
Lucid$treatlist[Lucid$sensitivelist==1] <- 1
Lucid$treatlist[Lucid$sensitivelist==2] <- 2
Lucid$treatlist[Lucid$sensitivelist==3] <- 3
Lucid$treatlist[Lucid$sensitivelist==4] <- 4

Lucid$list <- Lucid$controllist + Lucid$treatlist

##Create Treatment Variables
Lucid$treatment <- 0
Lucid$treatment[Lucid$treatment_variable==1] <- 1  

##Bootstrapped Lucid Results##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$list[which(temp$treatment==1)], na.rm=TRUE)
    B <- mean(temp$list[which(temp$treatment==0)], na.rm=TRUE)
    C <- mean(temp$direct, na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Lucid)


##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])


####################Appendix####################

##########Figure A.2##########
chem1 <- Qualtrics
data.list1 <- chem1[!(is.na(chem1$chem_group)),]

multi1 <- ictreg(chem_all ~ gender + age_cat + college + republican + democrat + income + white + black , data = data.list1, treat = "chem_group", J = 3, method = "ml", overdispersed = F, constrained = T)
summary(multi1)

gender.male <- gender.female <- data.list1
gender.male[, "gender"] <- 1
gender.female[, "gender"] <- 2

avg.pred.male.mle <- predict(multi1, newdata = gender.male, avg = T, se.fit = T, interval = "confidence")
avg.pred.female.mle <- predict(multi1, newdata = gender.female, avg = T, se.fit = T, interval = "confidence")

age.1829 <- age.3039 <- age.4049 <- age.5064 <- age.65 <- data.list1
age.1829[, "age_cat"] <- 1
age.3039[, "age_cat"] <- 2
age.4049[, "age_cat"] <- 3
age.5064[, "age_cat"] <- 4
age.65[, "age_cat"] <- 5

avg.pred.1829.mle <- predict(multi1, newdata = age.1829, avg = T, se.fit = T, interval = "confidence")
avg.pred.3039.mle <- predict(multi1, newdata = age.3039, avg = T, se.fit = T, interval = "confidence")
avg.pred.4049.mle <- predict(multi1, newdata = age.4049, avg = T, se.fit = T, interval = "confidence")
avg.pred.5064.mle <- predict(multi1, newdata = age.5064, avg = T, se.fit = T, interval = "confidence")
avg.pred.65.mle <- predict(multi1, newdata = age.65, avg = T, se.fit = T, interval = "confidence")

edu.nocoll <- edu.college <- data.list1
edu.nocoll[, "college"] <- 0
edu.college[, "college"] <- 1

avg.pred.nocoll.mle <- predict(multi1, newdata = edu.nocoll, avg = T, se.fit = T, interval = "confidence")
avg.pred.coll.mle <- predict(multi1, newdata = edu.college, avg = T, se.fit = T, interval = "confidence")

party.dem <- party.rep <- party.ind <- data.list1
party.dem[, "democrat"] <- 1
party.dem[, "republican"] <- 0
party.dem[, "independent"] <- 0
party.rep[, "democrat"] <- 0
party.rep[, "republican"] <- 1
party.rep[, "independent"] <- 0
party.ind[, "democrat"] <- 0
party.ind[, "republican"] <- 0
party.ind[, "independent"] <- 1

avg.pred.dem.mle <- predict(multi1, newdata = party.dem, avg = T, se.fit = T, interval = "confidence")
avg.pred.rep.mle <- predict(multi1, newdata = party.rep, avg = T, se.fit = T, interval = "confidence")
avg.pred.ind.mle <- predict(multi1, newdata = party.ind, avg = T, se.fit = T, interval = "confidence")

class.low <- class.mid <- class.upper <- data.list1
class.low[, "income"] <- 1
class.low[, "income"] <- 2
class.mid[, "income"] <- 3
class.mid[, "income"] <- 4
class.upper[, "income"] <- 5
class.upper[, "income"] <- 6
class.upper[, "income"] <- 7

avg.pred.low.mle <- predict(multi1, newdata = class.low, avg = T, se.fit = T, interval = "confidence")
avg.pred.mid.mle <- predict(multi1, newdata = class.mid, avg = T, se.fit = T, interval = "confidence")
avg.pred.upper.mle <- predict(multi1, newdata = class.upper, avg = T, se.fit = T, interval = "confidence")

race.other <- race.black <- race.white <- data.list1
race.white[, "white"] <- 1
race.black[, "black"] <- 1
race.other[, "white"] <- 0
race.other[, "black"] <- 0

avg.pred.white.mle <- predict(multi1, newdata = race.white, avg = T, se.fit = T, interval = "confidence")
avg.pred.black.mle <- predict(multi1, newdata = race.black, avg = T, se.fit = T, interval = "confidence")
avg.pred.other.mle <- predict(multi1, newdata = race.other, avg = T, se.fit = T, interval = "confidence")

summary(avg.pred.male.mle$fit)
summary(avg.pred.male.mle$se.fit)

summary(avg.pred.female.mle$fit)
summary(avg.pred.female.mle$se.fit)

summary(avg.pred.1829.mle$fit)
summary(avg.pred.1829.mle$se.fit)

summary(avg.pred.3039.mle$fit)
summary(avg.pred.3039.mle$se.fit)

summary(avg.pred.4049.mle$fit)
summary(avg.pred.4049.mle$se.fit)

summary(avg.pred.5064.mle$fit)
summary(avg.pred.5064.mle$se.fit)

summary(avg.pred.65.mle$fit)
summary(avg.pred.65.mle$se.fit)

summary(avg.pred.nocoll.mle$fit)
summary(avg.pred.nocoll.mle$se.fit)

summary(avg.pred.coll.mle$fit)
summary(avg.pred.coll.mle$se.fit)

summary(avg.pred.dem.mle$fit)
summary(avg.pred.dem.mle$se.fit)

summary(avg.pred.ind.mle$fit)
summary(avg.pred.ind.mle$se.fit)

summary(avg.pred.rep.mle$fit)
summary(avg.pred.rep.mle$se.fit)

summary(avg.pred.low.mle$fit)
summary(avg.pred.low.mle$se.fit)

summary(avg.pred.mid.mle$fit)
summary(avg.pred.mid.mle$se.fit)

summary(avg.pred.upper.mle$fit)
summary(avg.pred.upper.mle$se.fit)

summary(avg.pred.white.mle$fit)
summary(avg.pred.white.mle$se.fit)

summary(avg.pred.black.mle$fit)
summary(avg.pred.black.mle$se.fit)

summary(avg.pred.other.mle$fit)
summary(avg.pred.other.mle$se.fit)

coef.vec <- c(0.2973, 0.366, 0.2284, 0.2751, 0.3273, 0.384, 0.4439, 0.3452, 0.3069, 0.363, 0.3052, 0.3246, 0.3578, 0.313, 0.2519, 0.312, 0.3076, 0.4024)
se.vec <- c(0.1069, 0.07424, 0.08853, 0.07121, 0.06523, 0.08591, 0.1257, 0.08024, 0.1224, 0.1092, 0.09348, 0.1458, 0.08124, 0.07213, 0.1341, 0.07319, 0.2122, 0.2312)
var.names <- c("Male", "Female", "Age: 18-29", "Age: 30-39", "Age: 40-49", "Age: 50-64", "Age: 65+", "No College", "College", "Democrat", "Independent", "Republican", "Lower Class" , "Middle Class" , "Upper Class", "White", "Black", "Other Race")

y.axis <- length(var.names):1

layout(matrix(c(2,1),1,2), widths = c(1.5, 5))

par(mar=c(2,5,.5,1))

plot(coef.vec, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, 
     xlim = c(0,1), xaxs = "r", main = "")

segments(coef.vec-qnorm(.975)*se.vec, y.axis, coef.vec+qnorm(.975)*se.vec, y.axis, lwd =  1.5)

axis(1, at = seq(0,1,by=.1), labels =  seq(0,1,by=.1), tick = T,
     cex.axis = .8, mgp = c(2,.5,0))
axis(2, at = y.axis, label = var.names, las = 1, tick = T, 
     cex.axis = .8) 

par(mar=c(2,0,.5,0)) 

plot(seq(0,1,length=length(var.names)), y.axis, type = "n", axes = F,  xlab = "Estimated Proportion of Respondents Supporting Chemical Weapons Use", ylab = "")


##########Table A.3########## 
Qualtrics$ThirtytoThirtyNine <- 0
Qualtrics$ThirtytoThirtyNine[Qualtrics$age>=30 & Qualtrics$age<=39] <- 1

Qualtrics$FortytoFortyNine <- 0
Qualtrics$FortytoFortyNine[Qualtrics$age>=40 & Qualtrics$age<=49] <- 1

Qualtrics$FiftytoSixtyFour <- 0
Qualtrics$FiftytoSixtyFour[Qualtrics$age>=50 & Qualtrics$age<=64] <- 1

Qualtrics$SixtyFivePlus <- 0
Qualtrics$SixtyFivePlus[Qualtrics$age>=65] <- 1

Qualtrics$HundredtoHundredFifty <- 0
Qualtrics$HundredtoHundredFifty[Qualtrics$income==5] <- 1

Qualtrics$HundredFiftytoTwoHundred <- 0
Qualtrics$HundredFiftytoTwoHundred[Qualtrics$income==6] <- 1

Qualtrics$TwoHundredPlus <- 0
Qualtrics$TwoHundredPlus[Qualtrics$income==7] <- 1

Qualtrics$YesDirect <- NA
Qualtrics$YesDirect[Qualtrics$chem1==0] <- 0
Qualtrics$YesDirect[Qualtrics$chem1==1] <- 1

Qualtrics$NotSure <- NA
Qualtrics$NotSure[Qualtrics$chem1==0 | Qualtrics$chem1==1] <- 0
Qualtrics$NotSure[Qualtrics$chem1==4] <- 1

Qualtrics$YesOrNotSure <- NA
Qualtrics$YesOrNotSure[Qualtrics$chem1==0] <- 0
Qualtrics$YesOrNotSure[Qualtrics$chem1==1 | Qualtrics$chem1==4] <- 1

##Model 1##
model1 <- glm(YesDirect ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + college + republican + democrat +  TwentyFivetoFiftyK + FiftytoSeventyFiveK + SeventyFivetoHundredK + HundredtoHundredFifty + HundredFiftytoTwoHundred + TwoHundredPlus + White + Black, data=Qualtrics, family=binomial(link='logit'))

cov.m1 <- vcovHC(model1, type = "HC0")

std.err1 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model1)
  , "Robust SE" = std.err1
  , z = (coef(model1)/std.err1)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model1)/std.err1), lower.tail = FALSE)
  , LL = coef(model1) - q.val  * std.err1
  , UL = coef(model1) + q.val  * std.err1
)

r.est

##Model 2##
model2 <- glm(NotSure ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + college + republican + democrat +  TwentyFivetoFiftyK + FiftytoSeventyFiveK + SeventyFivetoHundredK + HundredtoHundredFifty + HundredFiftytoTwoHundred + TwoHundredPlus + White + Black, data=Qualtrics, family=binomial(link='logit'))

cov.m1 <- vcovHC(model2, type = "HC0")

std.err2 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model2)
  , "Robust SE" = std.err2
  , z = (coef(model2)/std.err2)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model2)/std.err2), lower.tail = FALSE)
  , LL = coef(model2) - q.val  * std.err2
  , UL = coef(model2) + q.val  * std.err2
)

r.est

##Model 3##
model3 <- glm(YesOrNotSure ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + college + republican + democrat +  TwentyFivetoFiftyK + FiftytoSeventyFiveK + SeventyFivetoHundredK + HundredtoHundredFifty + HundredFiftytoTwoHundred + TwoHundredPlus + White + Black, data=Qualtrics, family=binomial(link='logit'))

cov.m1 <- vcovHC(model3, type = "HC0")

std.err3 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model3)
  , "Robust SE" = std.err3
  , z = (coef(model3)/std.err3)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model3)/std.err3), lower.tail = FALSE)
  , LL = coef(model3) - q.val  * std.err3
  , UL = coef(model3) + q.val  * std.err3
)

r.est


##########Table A.5##########

##Dropping Attention Check Failures##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==1 & temp$attention==1)], temp$wtChem[which(temp$treatment==1 & temp$t_chem==1 & temp$attention==1)], na.rm=TRUE)
    B <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==1 & temp$attention==1)], temp$wtChem[which(temp$treatment==0 & temp$t_chem==1 & temp$attention==1)], na.rm=TRUE)
    C <- weighted.mean(temp$direct[which(temp$t_chem==1 & temp$attention==1)], temp$wtChem[which(temp$t_chem==1 & temp$attention==1)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.3 <- bootTreat2(Turk)

##Indirect Support##
round(c(full.3$Indirect, quantile(full.3$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,1] > 0))/length(full.3$boot[,1])

##Direct Support##
round(c(full.3$Direct, quantile(full.3$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,2] > 0))/length(full.3$boot[,2])

##Difference (Indirect-Direct)
round(c(full.3$Superficial, quantile(full.3$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,3] > 0))/length(full.3$boot[,3])

###Unweighted##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$list[which(temp$treatment==1 & temp$t_chem==1)], na.rm=TRUE)
    B <- mean(temp$list[which(temp$treatment==0 & temp$t_chem==1)], na.rm=TRUE)
    C <- mean(temp$direct[which(temp$t_chem==1)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Turk)

##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])


##########Table A.6##########

##Full Weighted Sample##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==0)], temp$wtNuke[which(temp$treatment==1 & temp$t_chem==0)], na.rm=TRUE)
    B <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==0)], temp$wtNuke[which(temp$treatment==0 & temp$t_chem==0)], na.rm=TRUE)
    C <- weighted.mean(temp$direct[which(temp$t_chem==0)], temp$wtNuke[which(temp$t_chem==0)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Turk)

##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])

##Dropping Attention Check Failures##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==0 & temp$attention==1)], temp$wtNuke[which(temp$treatment==1 & temp$t_chem==0 & temp$attention==1)], na.rm=TRUE)
    B <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==0 & temp$attention==1)], temp$wtNuke[which(temp$treatment==0 & temp$t_chem==0 & temp$attention==1)], na.rm=TRUE)
    C <- weighted.mean(temp$direct[which(temp$t_chem==0 & temp$attention==1)], temp$wtNuke[which(temp$t_chem==0 & temp$attention==1)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.3 <- bootTreat2(Turk)

##Indirect Support##
round(c(full.3$Indirect, quantile(full.3$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,1] > 0))/length(full.3$boot[,1])

##Direct Support##
round(c(full.3$Direct, quantile(full.3$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,2] > 0))/length(full.3$boot[,2])

##Difference (Indirect-Direct)
round(c(full.3$Superficial, quantile(full.3$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,3] > 0))/length(full.3$boot[,3])

##Unweighted##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$list[which(temp$treatment==1 & temp$t_chem==0)], na.rm=TRUE)
    B <- mean(temp$list[which(temp$treatment==0 & temp$t_chem==0)], na.rm=TRUE)
    C <- mean(temp$direct[which(temp$t_chem==0)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Turk)

##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])


##########Figure A.7##########
data.list2chem <- Turk[which(Turk$t_chem==1),]

multi2 <- ictreg(list ~ PeerCue + EliteCue + Female + age_cat + college + Republican + Democrat + income_cat + White + Black + RPI_Median, data = data.list2chem, treat = "treatment", weight="wtChem", J = 3, method = "ml", overdispersed = F, constrained = T)
summary(multi2)

cue.peer <- cue.elite <- cue.none <- data.list2chem
cue.none[, "NoCue"] <- 1
cue.peer[, "PeerCue"] <- 1
cue.elite[, "EliteCue"] <- 1

avg.pred.nocue.mle <- predict(multi2, newdata = cue.none, avg = T, se.fit = T, interval = "confidence")
avg.pred.peercue.mle <- predict(multi2, newdata = cue.peer, avg = T, se.fit = T, interval = "confidence")
avg.pred.elitecue.mle <- predict(multi2, newdata = cue.elite, avg = T, se.fit = T, interval = "confidence")

gender.male <- gender.female <- data.list2chem
gender.male[, "Female"] <- 0
gender.female[, "Female"] <- 1

avg.pred.male.mle <- predict(multi2, newdata = gender.male, avg = T, se.fit = T, interval = "confidence")
avg.pred.female.mle <- predict(multi2, newdata = gender.female, avg = T, se.fit = T, interval = "confidence")

age.1829 <- age.3039 <- age.4049 <- age.5064 <- age.65 <- data.list2chem
age.1829[, "age_cat"] <- 1
age.3039[, "age_cat"] <- 2
age.4049[, "age_cat"] <- 3
age.5064[, "age_cat"] <- 4
age.65[, "age_cat"] <- 5

avg.pred.1829.mle <- predict(multi2, newdata = age.1829, avg = T, se.fit = T, interval = "confidence")
avg.pred.3039.mle <- predict(multi2, newdata = age.3039, avg = T, se.fit = T, interval = "confidence")
avg.pred.4049.mle <- predict(multi2, newdata = age.4049, avg = T, se.fit = T, interval = "confidence")
avg.pred.5064.mle <- predict(multi2, newdata = age.5064, avg = T, se.fit = T, interval = "confidence")
avg.pred.65.mle <- predict(multi2, newdata = age.65, avg = T, se.fit = T, interval = "confidence")

edu.nocoll <- edu.college <- data.list2chem
edu.nocoll[, "college"] <- 0
edu.college[, "college"] <- 1

avg.pred.nocoll.mle <- predict(multi2, newdata = edu.nocoll, avg = T, se.fit = T, interval = "confidence")
avg.pred.coll.mle <- predict(multi2, newdata = edu.college, avg = T, se.fit = T, interval = "confidence")

party.dem <- party.rep <- party.ind <- data.list2chem
party.dem[, "Democrat"] <- 1
party.dem[, "Republican"] <- 0
party.dem[, "Independent"] <- 0
party.rep[, "Democrat"] <- 0
party.rep[, "Republican"] <- 1
party.rep[, "Independent"] <- 0
party.ind[, "Democrat"] <- 0
party.ind[, "Republican"] <- 0
party.ind[, "Independent"] <- 1

avg.pred.dem.mle <- predict(multi2, newdata = party.dem, avg = T, se.fit = T, interval = "confidence")
avg.pred.rep.mle <- predict(multi2, newdata = party.rep, avg = T, se.fit = T, interval = "confidence")
avg.pred.ind.mle <- predict(multi2, newdata = party.ind, avg = T, se.fit = T, interval = "confidence")

class.low <- class.mid <- class.upper <- data.list2chem
class.low[, "income_cat"] <- 1
class.low[, "income_cat"] <- 2
class.mid[, "income_cat"] <- 3
class.mid[, "income_cat"] <- 4
class.upper[, "income_cat"] <- 5
class.upper[, "income_cat"] <- 6

avg.pred.low.mle <- predict(multi2, newdata = class.low, avg = T, se.fit = T, interval = "confidence")
avg.pred.mid.mle <- predict(multi2, newdata = class.mid, avg = T, se.fit = T, interval = "confidence")
avg.pred.upper.mle <- predict(multi2, newdata = class.upper, avg = T, se.fit = T, interval = "confidence")

race.other <- race.black <- race.white <- data.list2chem
race.white[, "White"] <- 1
race.black[, "Black"] <- 1
race.other[, "Black"] <- 0
race.other[, "White"] <- 0

avg.pred.white.mle <- predict(multi2, newdata = race.white, avg = T, se.fit = T, interval = "confidence")
avg.pred.black.mle <- predict(multi2, newdata = race.black, avg = T, se.fit = T, interval = "confidence")
avg.pred.other.mle <- predict(multi2, newdata = race.other, avg = T, se.fit = T, interval = "confidence")

rpi.low <- rpi.high <- data.list2chem
rpi.high[, "RPI_Median"] <- 1
rpi.low[, "RPI_Median"] <- 0

avg.pred.rpilow.mle <- predict(multi2, newdata = rpi.low, avg = T, se.fit = T, interval = "confidence")
avg.pred.rpihigh.mle <- predict(multi2, newdata = rpi.high, avg = T, se.fit = T, interval = "confidence")

summary(avg.pred.nocue.mle$fit)
summary(avg.pred.nocue.mle$se.fit)

summary(avg.pred.peercue.mle$fit)
summary(avg.pred.peercue.mle$se.fit)

summary(avg.pred.elitecue.mle$fit)
summary(avg.pred.elitecue.mle$se.fit)

summary(avg.pred.male.mle$fit)
summary(avg.pred.male.mle$se.fit)

summary(avg.pred.female.mle$fit)
summary(avg.pred.female.mle$se.fit)

summary(avg.pred.1829.mle$fit)
summary(avg.pred.1829.mle$se.fit)

summary(avg.pred.3039.mle$fit)
summary(avg.pred.3039.mle$se.fit)

summary(avg.pred.4049.mle$fit)
summary(avg.pred.4049.mle$se.fit)

summary(avg.pred.5064.mle$fit)
summary(avg.pred.5064.mle$se.fit)

summary(avg.pred.65.mle$fit)
summary(avg.pred.65.mle$se.fit)

summary(avg.pred.nocoll.mle$fit)
summary(avg.pred.nocoll.mle$se.fit)

summary(avg.pred.coll.mle$fit)
summary(avg.pred.coll.mle$se.fit)

summary(avg.pred.dem.mle$fit)
summary(avg.pred.dem.mle$se.fit)

summary(avg.pred.ind.mle$fit)
summary(avg.pred.ind.mle$se.fit)

summary(avg.pred.rep.mle$fit)
summary(avg.pred.rep.mle$se.fit)

summary(avg.pred.low.mle$fit)
summary(avg.pred.low.mle$se.fit)

summary(avg.pred.mid.mle$fit)
summary(avg.pred.mid.mle$se.fit)

summary(avg.pred.upper.mle$fit)
summary(avg.pred.upper.mle$se.fit)

summary(avg.pred.white.mle$fit)
summary(avg.pred.white.mle$se.fit)

summary(avg.pred.black.mle$fit)
summary(avg.pred.black.mle$se.fit)

summary(avg.pred.other.mle$fit)
summary(avg.pred.other.mle$se.fit)

summary(avg.pred.rpilow.mle$fit)
summary(avg.pred.rpilow.mle$se.fit)

summary(avg.pred.rpihigh.mle$fit)
summary(avg.pred.rpihigh.mle$se.fit)

coef.vec <- c(0.07683, 0.06613, 0.1931, 0.1224, 0.03961, 0.04544, 0.06684, 0.1024, 0.1544, 0.2191, 0.1879, 0.02471, 0.02231, 0.06456, 0.2104, 0.04591, 0.119, 0.2531, 0.04846, 0.04115, 0.2514, 0.1289, 0.04819)
se.vec <- c(0.009716, 0.02085, 0.03572, 0.01586, 0.008111, 0.009, 0.01027, 0.01269, 0.01922, 0.03289, 0.02369, 0.007742, 0.009747, 0.01396, 0.02867, 0.009146, 0.01362, 0.03191, 0.01053, 0.01035, 0.03121, 0.02011, 0.007213)
var.names <- c("No Cue", "Peer Cue", "Elite Cue", "Male", "Female", "Age: 18-29", "Age: 30-39", "Age: 40-49", "Age: 50-64", "Age: 65+", "No College",  "College", "Democrat", "Independent", "Republican", "Low Class", "Middle Class", "Upper Class", "White", "Black", "Other Race", "Low RPI", "High RPI")

y.axis <- length(var.names):1

layout(matrix(c(2,1),1,2), widths = c(1.5, 5))

par(mar=c(2,5,.5,1))

plot(coef.vec, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, 
     xlim = c(0,1), xaxs = "r", main = "")

segments(coef.vec-qnorm(.975)*se.vec, y.axis, coef.vec+qnorm(.975)*se.vec, y.axis, lwd =  1.5)

axis(1, at = seq(0,1,by=.1), labels =  seq(0,1,by=.1), tick = T,
     cex.axis = .8, mgp = c(2,.5,0))
axis(2, at = y.axis, label = var.names, las = 1, tick = T, 
     cex.axis = .8) 

par(mar=c(2,0,.5,0)) 

plot(seq(0,1,length=length(var.names)), y.axis, type = "n", axes = F,  xlab = "Estimated Proportion of Respondents Supporting Chemical Weapons Use", ylab = "")


##########Figure A.8##########
data.list2nuke <- Turk[which(Turk$t_chem==0),]

multi3 <- ictreg(list ~ PeerCue + EliteCue + Female + age_cat + college + Republican + Democrat + income_cat + White + Black + RPI_Median, data = data.list2nuke, treat = "treatment", weight="wtNuke", J = 3, method = "ml", overdispersed = F, constrained = T)
summary(multi3)

cue.peer <- cue.elite <- cue.none <- data.list2nuke
cue.none[, "NoCue"] <- 1
cue.peer[, "PeerCue"] <- 1
cue.elite[, "EliteCue"] <- 1

avg.pred.nocue.mle <- predict(multi3, newdata = cue.none, avg = T, se.fit = T, interval = "confidence")
avg.pred.peercue.mle <- predict(multi3, newdata = cue.peer, avg = T, se.fit = T, interval = "confidence")
avg.pred.elitecue.mle <- predict(multi3, newdata = cue.elite, avg = T, se.fit = T, interval = "confidence")

gender.male <- gender.female <- data.list2nuke
gender.male[, "Female"] <- 0
gender.female[, "Female"] <- 1

avg.pred.male.mle <- predict(multi3, newdata = gender.male, avg = T, se.fit = T, interval = "confidence")
avg.pred.female.mle <- predict(multi3, newdata = gender.female, avg = T, se.fit = T, interval = "confidence")

age.1829 <- age.3039 <- age.4049 <- age.5064 <- age.65 <- data.list2nuke
age.1829[, "age_cat"] <- 1
age.3039[, "age_cat"] <- 2
age.4049[, "age_cat"] <- 3
age.5064[, "age_cat"] <- 4
age.65[, "age_cat"] <- 5

avg.pred.1829.mle <- predict(multi3, newdata = age.1829, avg = T, se.fit = T, interval = "confidence")
avg.pred.3039.mle <- predict(multi3, newdata = age.3039, avg = T, se.fit = T, interval = "confidence")
avg.pred.4049.mle <- predict(multi3, newdata = age.4049, avg = T, se.fit = T, interval = "confidence")
avg.pred.5064.mle <- predict(multi3, newdata = age.5064, avg = T, se.fit = T, interval = "confidence")
avg.pred.65.mle <- predict(multi3, newdata = age.65, avg = T, se.fit = T, interval = "confidence")

edu.nocoll <- edu.college <- data.list2nuke
edu.nocoll[, "college"] <- 0
edu.college[, "college"] <- 1

avg.pred.nocoll.mle <- predict(multi3, newdata = edu.nocoll, avg = T, se.fit = T, interval = "confidence")
avg.pred.coll.mle <- predict(multi3, newdata = edu.college, avg = T, se.fit = T, interval = "confidence")

party.dem <- party.rep <- party.ind <- data.list2nuke
party.dem[, "Democrat"] <- 1
party.dem[, "Republican"] <- 0
party.dem[, "Independent"] <- 0
party.rep[, "Democrat"] <- 0
party.rep[, "Republican"] <- 1
party.rep[, "Independent"] <- 0
party.ind[, "Democrat"] <- 0
party.ind[, "Republican"] <- 0
party.ind[, "Independent"] <- 1

avg.pred.dem.mle <- predict(multi3, newdata = party.dem, avg = T, se.fit = T, interval = "confidence")
avg.pred.rep.mle <- predict(multi3, newdata = party.rep, avg = T, se.fit = T, interval = "confidence")
avg.pred.ind.mle <- predict(multi3, newdata = party.ind, avg = T, se.fit = T, interval = "confidence")

class.low <- class.mid <- class.upper <- data.list2nuke
class.low[, "income_cat"] <- 1
class.low[, "income_cat"] <- 2
class.mid[, "income_cat"] <- 3
class.mid[, "income_cat"] <- 4
class.upper[, "income_cat"] <- 5
class.upper[, "income_cat"] <- 6

avg.pred.low.mle <- predict(multi3, newdata = class.low, avg = T, se.fit = T, interval = "confidence")
avg.pred.mid.mle <- predict(multi3, newdata = class.mid, avg = T, se.fit = T, interval = "confidence")
avg.pred.upper.mle <- predict(multi3, newdata = class.upper, avg = T, se.fit = T, interval = "confidence")

race.other <- race.black <- race.white <- data.list2nuke
race.white[, "White"] <- 1
race.black[, "Black"] <- 1
race.other[, "Black"] <- 0
race.other[, "White"] <- 0

avg.pred.white.mle <- predict(multi3, newdata = race.white, avg = T, se.fit = T, interval = "confidence")
avg.pred.black.mle <- predict(multi3, newdata = race.black, avg = T, se.fit = T, interval = "confidence")
avg.pred.other.mle <- predict(multi3, newdata = race.other, avg = T, se.fit = T, interval = "confidence")

rpi.low <- rpi.high <- data.list2nuke
rpi.high[, "RPI_Median"] <- 1
rpi.low[, "RPI_Median"] <- 0

avg.pred.rpilow.mle <- predict(multi3, newdata = rpi.low, avg = T, se.fit = T, interval = "confidence")
avg.pred.rpihigh.mle <- predict(multi3, newdata = rpi.high, avg = T, se.fit = T, interval = "confidence")

summary(avg.pred.nocue.mle$fit)
summary(avg.pred.nocue.mle$se.fit)

summary(avg.pred.peercue.mle$fit)
summary(avg.pred.peercue.mle$se.fit)

summary(avg.pred.elitecue.mle$fit)
summary(avg.pred.elitecue.mle$se.fit)

summary(avg.pred.male.mle$fit)
summary(avg.pred.male.mle$se.fit)

summary(avg.pred.female.mle$fit)
summary(avg.pred.female.mle$se.fit)

summary(avg.pred.1829.mle$fit)
summary(avg.pred.1829.mle$se.fit)

summary(avg.pred.3039.mle$fit)
summary(avg.pred.3039.mle$se.fit)

summary(avg.pred.4049.mle$fit)
summary(avg.pred.4049.mle$se.fit)

summary(avg.pred.5064.mle$fit)
summary(avg.pred.5064.mle$se.fit)

summary(avg.pred.65.mle$fit)
summary(avg.pred.65.mle$se.fit)

summary(avg.pred.nocoll.mle$fit)
summary(avg.pred.nocoll.mle$se.fit)

summary(avg.pred.coll.mle$fit)
summary(avg.pred.coll.mle$se.fit)

summary(avg.pred.dem.mle$fit)
summary(avg.pred.dem.mle$se.fit)

summary(avg.pred.ind.mle$fit)
summary(avg.pred.ind.mle$se.fit)

summary(avg.pred.rep.mle$fit)
summary(avg.pred.rep.mle$se.fit)

summary(avg.pred.low.mle$fit)
summary(avg.pred.low.mle$se.fit)

summary(avg.pred.mid.mle$fit)
summary(avg.pred.mid.mle$se.fit)

summary(avg.pred.upper.mle$fit)
summary(avg.pred.upper.mle$se.fit)

summary(avg.pred.white.mle$fit)
summary(avg.pred.white.mle$se.fit)

summary(avg.pred.black.mle$fit)
summary(avg.pred.black.mle$se.fit)

summary(avg.pred.other.mle$fit)
summary(avg.pred.other.mle$se.fit)

summary(avg.pred.rpilow.mle$fit)
summary(avg.pred.rpilow.mle$se.fit)

summary(avg.pred.rpihigh.mle$fit)
summary(avg.pred.rpihigh.mle$se.fit)

coef.vec <- c(0.5528, 0.4254, 0.361, 0.4555, 0.6277, 0.5456, 0.5513, 0.557, 0.5627, 0.5684, 0.1807, 0.8803, 0.6777, 0.4967, 0.3776, 0.5589, 0.5453, 0.5316, 0.6236, 0.7326, 0.3117, 0.5909, 0.5094)
se.vec <- c(0.0375, 0.09081, 0.09136, 0.0476, 0.03959, 0.05806, 0.04095, 0.03223, 0.03782, 0.05318, 0.05467, 0.0581, 0.05869, 0.0789, 0.06644, 0.04396, 0.03782, 0.06095, 0.0463, 0.06481, 0.05401, 0.05069, 0.04863)
var.names <- c("No Cue", "Peer Cue", "Elite Cue", "Male", "Female", "Age: 18-29", "Age: 30-39", "Age: 40-49", "Age: 50-64", "Age: 65+", "No College",  "College", "Democrat", "Independent", "Republican", "Low Class", "Middle Class", "Upper Class", "White", "Black", "Other Race", "Low RPI", "High RPI")

y.axis <- length(var.names):1

layout(matrix(c(2,1),1,2), widths = c(1.5, 5))

par(mar=c(2,5,.5,1))

plot(coef.vec, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, 
     xlim = c(0,1), xaxs = "r", main = "")

segments(coef.vec-qnorm(.975)*se.vec, y.axis, coef.vec+qnorm(.975)*se.vec, y.axis, lwd =  1.5)

axis(1, at = seq(0,1,by=.1), labels =  seq(0,1,by=.1), tick = T,
     cex.axis = .8, mgp = c(2,.5,0))
axis(2, at = y.axis, label = var.names, las = 1, tick = T, 
     cex.axis = .8) 

par(mar=c(2,0,.5,0)) 

plot(seq(0,1,length=length(var.names)), y.axis, type = "n", axes = F,  xlab = "Estimated Proportion of Respondents Supporting Chemical Weapons Use", ylab = "")


##########Figure A.9########## 

Open_Friends_Chem <- subset(Open_Friends, t_chem!=0)
Open_Friends_Chem$text <- Open_Friends_Chem$open_friends
Open_Friends_Chem$text <- as.character(Open_Friends_Chem$text)
Open_Friends_Chem$Negative_label <- "Friends/Family Support Chemical Weapons Use"
Open_Friends_Chem$Negative_label[Open_Friends_Chem$Negative==1] <- "Friends/Family Oppose Chemical Weapons Use"
corpus <- corpus(Open_Friends_Chem)
dfm <- dfm(corpus_subset(corpus, Negative_label %in% c("Friends/Family Oppose Chemical Weapons Use", "Friends/Family Support Chemical Weapons Use")), remove=c((stopwords("english")), "weapons"),
           remove_punct=TRUE, remove_numbers=TRUE, remove_symbols=TRUE, groups="Negative_label") %>% dfm_trim(min_termfreq = 3)
textplot_wordcloud(dfm, comparison=TRUE, max_words=300, color=c("blue", "red"))

##########Table A.10##########

#####Direct Support#####
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$direct[which(temp$t_chem==1 & temp$t_cue==0)], temp$wtChemNoCue[which(temp$t_chem==1 & temp$t_cue==0)], na.rm=TRUE)
    B <- weighted.mean(temp$direct[which(temp$t_chem==1 & temp$t_cue==1)], temp$wtChemPeerCue[which(temp$t_chem==1 & temp$t_cue==1)], na.rm=TRUE)
    C <- weighted.mean(temp$direct[which(temp$t_chem==1 & temp$t_cue==2)], temp$wtChemEliteCue[which(temp$t_chem==1 & temp$t_cue==2)], na.rm=TRUE)
    bootResults[i,1] <- (A) 
    bootResults[i,2] <- (B)
    bootResults[i,3] <- (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, NoCue=mean(bootResults[,1], na.rm=TRUE), 
              PeerCue=mean(bootResults[,2], na.rm=TRUE), 
              EliteCue=mean(bootResults[,3], na.rm=TRUE)))
}

full.3 <- bootTreat2(Turk)

##No Cue##
round(c(full.3$NoCue, quantile(full.3$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,1] > 0))/length(full.3$boot[,1])

##Peer Cue##
round(c(full.3$PeerCue, quantile(full.3$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,2] > 0))/length(full.3$boot[,2])

##Elite Cue##
round(c(full.3$EliteCue, quantile(full.3$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,3] > 0))/length(full.3$boot[,3])

##P-Value of Peer vs. No Cue##
1-length(which(full.3$boot[,2] > full.3$boot[,1]))/length(full.3$boot[,2])

##P-Value of Elite vs. No Cue##
1-length(which(full.3$boot[,3] > full.3$boot[,1]))/length(full.3$boot[,2])

#####List Support#####
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==1 & temp$t_cue==0)], temp$wtChemNoCue[which(temp$treatment==1 & temp$t_chem==1 & temp$t_cue==0)], na.rm=TRUE)
    B <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==1 & temp$t_cue==0)], temp$wtChemNoCue[which(temp$treatment==0 & temp$t_chem==1 & temp$t_cue==0)], na.rm=TRUE)
    C <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==1 & temp$t_cue==1)], temp$wtChemPeerCue[which(temp$treatment==1 & temp$t_chem==1 & temp$t_cue==1)], na.rm=TRUE)
    D <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==1 & temp$t_cue==1)], temp$wtChemPeerCue[which(temp$treatment==0 & temp$t_chem==1 & temp$t_cue==1)], na.rm=TRUE)
    E <- weighted.mean(temp$list[which(temp$treatment==1 & temp$t_chem==1 & temp$t_cue==2)], temp$wtChemEliteCue[which(temp$treatment==1 & temp$t_chem==1 & temp$t_cue==2)], na.rm=TRUE)
    G <- weighted.mean(temp$list[which(temp$treatment==0 & temp$t_chem==1 & temp$t_cue==2)], temp$wtChemEliteCue[which(temp$treatment==0 & temp$t_chem==1 & temp$t_cue==2)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C-D) 
    bootResults[i,3] <- (E-G) 
    drop(list())
  }
  return(list(model=label, boot=bootResults, NoCue=mean(bootResults[,1], na.rm=TRUE), 
              PeerCue=mean(bootResults[,2], na.rm=TRUE), 
              EliteCue=mean(bootResults[,3], na.rm=TRUE)))
}

full.3 <- bootTreat2(Turk)

##No Cue##
round(c(full.3$NoCue, quantile(full.3$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,1] > 0))/length(full.3$boot[,1])

##Peer Cue##
round(c(full.3$PeerCue, quantile(full.3$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,2] > 0))/length(full.3$boot[,2])

##Elite Cue##
round(c(full.3$EliteCue, quantile(full.3$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.3$boot[,3] > 0))/length(full.3$boot[,3])

##P-Value of Peer vs. No Cue##
1-length(which(full.3$boot[,2] > full.3$boot[,1]))/length(full.3$boot[,2])

##P-Value of Elite vs. No Cue##
1-length(which(full.3$boot[,3] > full.3$boot[,1]))/length(full.3$boot[,2])


##########Table A.11########## 
Turk$ThirtytoThirtyNine <- 0
Turk$ThirtytoThirtyNine[Turk$age>=30 & Turk$age<=39] <- 1

Turk$FortytoFortyNine <- 0
Turk$FortytoFortyNine[Turk$age>=40 & Turk$age<=49] <- 1

Turk$FiftytoSixtyFour <- 0
Turk$FiftytoSixtyFour[Turk$age>=50 & Turk$age<=64] <- 1

Turk$SixtyFivePlus <- 0
Turk$SixtyFivePlus[Turk$age>=65] <- 1

Turk$CollegeEducated <- 0
Turk$CollegeEducated[Turk$Education>=5] <- 1

Turk$notsure <- NA
Turk$notsure[Turk$chem_direct=="Yes"] <- 0
Turk$notsure[Turk$chem_direct=="No"] <- 0
Turk$notsure[Turk$chem_direct=="Not Sure"] <- 1

Turk$directnotsure <- NA
Turk$directnotsure[Turk$chem_direct=="No"] <- 0
Turk$directnotsure[Turk$chem_direct=="Not Sure" | Turk$chem_direct=="Yes"] <- 1

##Model 1##
model1 <- glm(direct ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + CollegeEducated + Republican + Democrat +  TwentytoFourtyFiveK + FourtyFivetoSeventyK + SeventytoHundredK + HundredtoTwoHundredK + TwoHundredK + White + Black + RPI_Median, weights=wtChem, data=subset(Turk, t_chem==1), family=binomial(link='logit'))

cov.m1 <- vcovHC(model1, type = "HC0")

std.err1 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model1)
  , "Robust SE" = std.err1
  , z = (coef(model1)/std.err1)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model1)/std.err1), lower.tail = FALSE)
  , LL = coef(model1) - q.val  * std.err1
  , UL = coef(model1) + q.val  * std.err1
)

r.est

##Model 2##
model2 <- glm(notsure ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + CollegeEducated + Republican + Democrat +  TwentytoFourtyFiveK + FourtyFivetoSeventyK + SeventytoHundredK + HundredtoTwoHundredK + TwoHundredK + White + Black + RPI_Median, weights=wtChem, data=subset(Turk, t_chem==1), family=binomial(link='logit'))

cov.m1 <- vcovHC(model2, type = "HC0")

std.err2 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model2)
  , "Robust SE" = std.err2
  , z = (coef(model2)/std.err2)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model2)/std.err2), lower.tail = FALSE)
  , LL = coef(model2) - q.val  * std.err2
  , UL = coef(model2) + q.val  * std.err2
)

r.est

##Model 3##
model3 <- glm(directnotsure ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + CollegeEducated + Republican + Democrat +  TwentytoFourtyFiveK + FourtyFivetoSeventyK + SeventytoHundredK + HundredtoTwoHundredK + TwoHundredK + White + Black + RPI_Median, weights=wtChem, data=subset(Turk, t_chem==1), family=binomial(link='logit'))

cov.m1 <- vcovHC(model3, type = "HC0")

std.err3 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model3)
  , "Robust SE" = std.err3
  , z = (coef(model3)/std.err3)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model3)/std.err3), lower.tail = FALSE)
  , LL = coef(model3) - q.val  * std.err3
  , UL = coef(model3) + q.val  * std.err3
)

r.est


##Model 4##
model4 <- glm(direct ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + CollegeEducated + Republican + Democrat +  TwentytoFourtyFiveK + FourtyFivetoSeventyK + SeventytoHundredK + HundredtoTwoHundredK + TwoHundredK + White + Black + RPI_Median, weights=wtNuke, data=subset(Turk, t_chem==0), family=binomial(link='logit'))

cov.m1 <- vcovHC(model4, type = "HC0")

std.err4 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model4)
  , "Robust SE" = std.err4
  , z = (coef(model4)/std.err4)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model4)/std.err4), lower.tail = FALSE)
  , LL = coef(model4) - q.val  * std.err4
  , UL = coef(model4) + q.val  * std.err4
)

r.est

##Model 5##
model5 <- glm(notsure ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + CollegeEducated + Republican + Democrat +  TwentytoFourtyFiveK + FourtyFivetoSeventyK + SeventytoHundredK + HundredtoTwoHundredK + TwoHundredK + White + Black + RPI_Median, weights=wtNuke, data=subset(Turk, t_chem==0), family=binomial(link='logit'))

cov.m1 <- vcovHC(model5, type = "HC0")

std.err5 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model5)
  , "Robust SE" = std.err5
  , z = (coef(model5)/std.err5)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model5)/std.err5), lower.tail = FALSE)
  , LL = coef(model5) - q.val  * std.err5
  , UL = coef(model5) + q.val  * std.err5
)

r.est

##Model 6##
model6 <- glm(directnotsure ~ Female + ThirtytoThirtyNine + FortytoFortyNine + FiftytoSixtyFour + SixtyFivePlus + CollegeEducated + Republican + Democrat +  TwentytoFourtyFiveK + FourtyFivetoSeventyK + SeventytoHundredK + HundredtoTwoHundredK + TwoHundredK + White + Black + RPI_Median, weights=wtNuke, data=subset(Turk, t_chem==0), family=binomial(link='logit'))

cov.m1 <- vcovHC(model6, type = "HC0")

std.err6 <- sqrt(diag(cov.m1))

q.val <- qnorm(0.975)

r.est <- cbind(
  Estimate = coef(model6)
  , "Robust SE" = std.err6
  , z = (coef(model6)/std.err6)
  , "Pr(>|z|) "= 2 * pnorm(abs(coef(model6)/std.err6), lower.tail = FALSE)
  , LL = coef(model6) - q.val  * std.err6
  , UL = coef(model6) + q.val  * std.err6
)

r.est


##########Table A.12########## 

data.list1 <- Qualtrics[!(is.na(Qualtrics$chem_group)),]
data.list2chem <- Turk[which(Turk$t_chem==1),]
data.list2nuke <- Turk[which(Turk$t_chem==0),]

test.value.study1 <- ict.test(data.list1$chem_all, data.list1$chem_group, J = 3, alpha = 0.05,  gms = T)
print(test.value.study1)

test.value.study2chem <- ict.test(data.list2chem$list, data.list2chem$treatment, J = 3, alpha = 0.05,  gms = T)
print(test.value.study2chem)

test.value.study2nuke <- ict.test(data.list2nuke$list, data.list2nuke$treatment, J = 3, alpha = 0.05,  gms = T)
print(test.value.study2nuke)


##########Table A.17########## 

###Unweighted##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$list[which(temp$treatment==1)], na.rm=TRUE)
    B <- mean(temp$list[which(temp$treatment==0)], na.rm=TRUE)
    C <- mean(temp$direct, na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Lucid)

##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])


##Dropping Attention Check Failures##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$list[which(temp$treatment==1 & temp$nuclearattention==1)], na.rm=TRUE)
    B <- mean(temp$list[which(temp$treatment==0 & temp$nuclearattention==1)], na.rm=TRUE)
    C <- mean(temp$direct[which(temp$nuclearattention==1)], na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Lucid)

##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])


##Weighted (All Respondents)##
set.seed(43215)
B <- 5000 

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- weighted.mean(temp$list[which(temp$treatment==1)], temp$wt[which(temp$treatment==1)], na.rm=TRUE)
    B <- weighted.mean(temp$list[which(temp$treatment==0)], temp$wt[which(temp$treatment==0)], na.rm=TRUE)
    C <- weighted.mean(temp$direct, temp$wt, na.rm=TRUE)
    bootResults[i,1] <- (A-B) 
    bootResults[i,2] <- (C)
    bootResults[i,3] <- (A-B) - (C)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Indirect=mean(bootResults[,1], na.rm=TRUE), 
              Direct=mean(bootResults[,2], na.rm=TRUE), 
              Superficial=mean(bootResults[,3], na.rm=TRUE)))
}

full.2 <- bootTreat2(Lucid)

##Indirect Support##
round(c(full.2$Indirect, quantile(full.2$boot[,1], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,1] > 0))/length(full.2$boot[,1])

##Direct Support##
round(c(full.2$Direct, quantile(full.2$boot[,2], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,2] > 0))/length(full.2$boot[,2])

##Difference (Indirect-Direct)
round(c(full.2$Superficial, quantile(full.2$boot[,3], c(0.025, 0.975))), digits=3)
#P-Value#
1-length(which(full.2$boot[,3] > 0))/length(full.2$boot[,3])








